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Abstract 


Non-charge  conserving  current  collection  algorithms  for  electromagnetic  PIC  plasma  simulations  may 
cause  errors  in  Gauss’s  law.  These  errors  arise  from  violations  of  the  charge  continuity  equation,  V  •  J  = 
—dp/dt^  which  in  turn  cause  errors  in  the  irrotational  part  of  E. 

Two  techniques  for  reducing  these  errors  are  examined  and  compared:  a  modified  Marder  correction 
which  corrects  electric  fields  locally  and  primarily  affects  short  wavelengths,  and  a  Boris  divergence 
correction,  which  solves  Poisson’s  equation  to  correct  the  electric  fields  so  that  Gauss’s  law  is  enforced 
globally.  The  effect  of  each  method  on  the  spectrum  of  the  error  is  examined.  Computational  efficiency 
and  accuracy  of  the  two  techniques  are  compared:  neither  method  is  clearly  superior. 

Cases  examined  include  corrections  in  electromagnetic  relativistic  beam  simulations,  and  a  hot  ther¬ 
mal  plasma.  In  addition,  the  spectral  comparison  provides  insight  into  the  behavior  of  the  schemes 
applied. 
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1  Introduction 


In  electromagnetic  particle-in-cell  (PIC)  simulations,  a  number  of  approximations  have  to  be  made  in  order 
to  model  the  physics  on  a  digital  computer  in  a  reasonable  time.  These  approximations  introduce  errors  in 
the  simulation  which  could  affect  the  physics. 

The  particular  approximations  which  are  commonly  made  include:  discretizing  space  for  the  purposes  of 
solution  of  Maxwell’s  curl  equations  (V  x  H  =  J  -f-  9D/5i,  V  x  E  =  — 9B/9t)  on  a  mathematical  mesh;  the 
use  of  many  fewer  particles  than  would  exist  in  a  real  system,  although  the  charge  to  mass  ratio  is  the  same; 
weighting  of  current  to  the  mathematical  mesh  from  these  discrete  particles  in  order  to  form  the  source  term 
J  for  solving  the  curl  equations;  and  the  weighting  of  the  fields  from  the  mathematical  mesh  to  the  particles 
to  calculate  their  motion.  These  approximations  are  typical  of  classical  PIC  models,  such  as  in  [1]  and  [2]. 

All  of  these  approximations  lead  to  error:  the  goal  is  to  keep  the  error  small  enough  so  that  the  simulation 
is  modeling  the  desired  physics  appropriately.  The  focus  of  this  work  is  to  study  the  accumulation  of  errors 
resulting  from  mapping  the  current  of  discrete  charged  particles  onto  a  numerical  mesh,  and  methods  of 
reducing  these  errors.  The  errors  arise  when  the  current  and  charge  density  mapped  to  the  mesh  violate  the 
continuity  equation.  A  number  of  authors  have  noted  this  problem  [1,  3,  4]. 

Non-charge  conserving  current  weightings,  such  as  a  linear  current  weighting  used  in  combination  with  a 
linear  charge  weighting,  violate  the  continuity  equation,  V- J  =  —dpldt,  in  turn  causing  errors  in  Gauss’s  law, 
V  •  D  =  p.  However,  these  linear  current  weightings  reduce  noise  in  comparison  with  the  charge  conserving 
method  described  by  Villasenor  and  Buneman  [5].  When  noise  is  an  issue,  one  may  choose  to  use  the  higher 
order  weighting  in  conjunction  with  a  correction  scheme  in  order  to  reduce  the  violation  of  Gauss’s  law. 
However,  higher  order  weightings  are  more  computationally  intensive  and  complicate  boundary  conditions. 

Presented  here  are  some  methods  which  can  be  used  to  reduce  cumulative  errors  in  Gauss’s  law  due  to 
slightly  non-ch2irge-conserving  current  weightings  from  the  particles  to  the  mesh,  and  the  results  observed 
when  these  methods  are  applied  to  a  variety  of  problems. 

Section  1.1  defines  the  corrections  used,  and  Section  1.2  describes  the  simulation  and  analysis  methods. 
Section  2.1  gives  results  for  removal  of  a  delta  function  source  test  case,  Section  2.2  for  single- wavelength 
errors.  In  Sections  2.3  and  2.4  results  for  relativistic  beams  are  presented.  Section  2.5  shows  results  for  a 
very  hot  thermal  plasma.  Conclusions  are  given  in  Section  3. 

1.1  Corrections  Used 

1.1.1  Boris  Correction 

The  Boris-type  correction  (modified  slightly  from  [1])  is  given  by: 


^corrected  —  E  “  WScj) 


(1) 
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such  that 


V  •  cEcorrected  —  P  (2) 

where  'Ecorrected  is  an  electric  field  which  obeys  Gauss’s  law,  and  V5<f)  is  the  correction  to  E,  the  uncorrected 
field.  Combining  Eqs.  (1)  and  (2),  we  obtain: 

V  •  (eE  -  =  p  (3) 

so  we  require  a  solution  for  5^,  with  E  and  p  given,  to  the  elliptic  equation: 

V  •  =  V  •  eE  -  p.  (4) 

This  elliptic  equation  can  be  solved  using  a  number  of  direct  and  iterative  schemes  with  varying  degrees 
of  efficiency. 


1.1.2  Langdon-Marder  Correction 


The  Marder-type  scheme  used  is  the  one  proposed  by  A.B.  Langdon  [3].  He  proposed  essentially  the  following 
scheme,  which  has  been  slightly  modified  to  allow  for  inhomogeneous  permittivity,  e: 

^Kt\ecte<i  =  +  AtV[d(V  •  (5) 

This  is  a  modification  of  the  Marder  scheme,  which  is  less  implicit: 


+  AtV[d(V  •  €E”  -  p”)] 


(6) 


where  d  is  the  diffusion  parameter,  which  should  satisfy  (in  two  dimensions' 


d< 


1  ,  Ax^Ay^  ^  _  j 


“  2 At  ^  Ax^  +  Ay2  ^ 


(7) 


to  ensure  stability  of  the  method  [3].  The  superscripts  n,  n  +  1  indicate  simulation  timestep:  Langdon’s 
proposed  modification  is  simply  to  use  the  E  and  p  at  the  current  time  rather  than  from  the  previous 
timestep. 


1.1.3  Irrotationality  of  the  correction 

Both  the  Boris  and  the  Langdon-Marder  corrections  arise  from  the  gradient  of  a  scalar  field,  V(5<^,  which  can 
be  seen  from  Equations  (4)  and  (5).  This,  combined  with  the  vector  identity  V  x  V5(^  =  0,  demonstrates 
that  both  corrections  have  no  effect  on  the  solenoidal  part  of  the  field.  This  property  is  important,  since  only 
the  solenoidal  part  of  the  electric  field  is  used  in  advancing  the  magnetic  field  in  time  via  V  x  E  =  — 
hence  both  methods  only  operate  on  the  irrotational  part  of  E. 
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1.2  Models  and  Implementation 


The  methods  were  implemented  and  tested  within  the  framework  of  the  XOOPIC  two-dimensional  electro¬ 
magnetic  PIC  code  [6].  For  this  study,  standard  PIC  methods  [1]  were  used  for  advancing  particles  and 
fields.  In  the  interest  of  clarity  and  convenience  of  Fourier  analysis,  this  investigation  is  performed  on  a 
uniform  Cartesian  mesh.  The  testing  space  was  discretized  into  either  a  16x16,  32x32,  or  64x64  square  mesh 
bounded  by  conductors. 

Figure  1  illustrates  the  locations  of  the  field  components  on  the  mesh  (which  is  a  standard  Yee  mesh  [7]), 
Electric  field  E  and  current  I  are  defined  at  the  midpoints  of  lines  connecting  the  nodes  of  the  mesh,  while 
p  and  <p  sxe  defined  on  the  nodes. 

The  Boris  correction  was  implemented  by  using  the  Dynamic  Alternating  Direction  Implicit  (DADI)  [8] 
method  to  solve  Eq.  (4),  which  is  discretized  using  the  standard  5-point  stencil  for  <j): 

~~  —  2^t,j  /g\ 

The  DADI  method  is  an  iterative  scheme  which  makes  use  of  a  “fictitious  timestep”  to  attempt  to  speed  up 
convergence  to  a  correct  solution.  Successive  iterations  of  DADI  will  choose  successive  “fictitious  timesteps” 
using  a  heuristic  rule.  If  a  poor  choice  of  initial  “fictitious  timestep”  is  made,  the  problem  will  take  longer 
to  converge  to  the  solution. 

The  Langdon-Marder  correction  was  implemented  by  casting  Eq.  (5)  into  finite  difference  form  as  in  the 
following  example  for  one  component  of  E: 


x,j+l/2,fc  corrected 


where  Err  is  the  divergence  error,  calculated  by: 


(9) 


—  ^y{^^y,j,k+l/2  ^  ^^yJ,k-l/2)  +  ^^(eEy ^^yJ,k^l/2)  “  Pjyk  (10) 


The  wavelength  components  of  the  error  are  analyzed  using  a  one  dimensional  Fast  Fourier  Transform 
along  each  mesh-line,  producing  kx  vs.  y,  and  ky  vs.  x.  Only  a  finite  set  of  wavelengths  are  observable: 
those  for  which  0  <  fegjAx  <  tt. 


5 


Figure  1:  Yee  mesh  for  current  and  electric  field. 


2  Results 


2.1  Delta  Function  Source  Test  Case 

The  methods  were  tested  first  on  a  delta  function  source,  by  setting  E*“®(x,2/)  =  0  and  setting  p  —  5{x  — 
—  Vo)-  A  delta  function  has  the  property  of  having  components  at  every  wavelength.  This  allows 
examination  of  the  behavior  of  each  correction  on  all  wavelengths  simultaneously. 

For  this  case,  the  bounded,  square  system  is  loaded  in  the  center  with  the  closest  approximation  to  a 
delta  function  source  possible  on  a  discrete  mesh:  a  nonzero  charge  density  at  only  one  mesh  point.  The 
walls  are  grounded  ideal  conductors  on  which  potentials  satisfy  a  Dirichlet  condition.  All  corrections  should 
relax  to  the  Poisson  solution  in  the  system  for  a  line  charge,  or  equivalently,  reduce  the  error  in  Gauss’s  law 
to  zero. 
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2.1.1  Boris  Correction  on  the  Delta  Function  Source 


Applying  the  Boris  correction  to  a  delta  function  source  shows  the  wavenumber  dependence  of  the  error  as 
a  function  of  iteration.  The  scales  of  Figure  2  are  logarithmic  over  9  orders  of  magnitude,  so  a  noticeable 
drop  in  amplitude  is  actually  fairly  large. 

The  error  f{k)  starts  at  a  (normalized)  value  of  1  for  all  wavelengths.  The  DADI-based  Boris  correction 
-requires  a  few  iterations  to  “settle  in” ,  and  find  the  correct  internal  parameter  for  the  most  rapid  convergence 
to  the  proper  solution.  The  parameter  is  a  “fictitious  timestep” ;  it  is  adjusted  at  each  iteration  by  the  DADI 
algorithm. 

The  first  try  at  reducing  the  error  shows  that  a  non-optimal  initial  “fictitious  timestep”  was  chosen.  Little 
removal  of  the  error  occurred  for  the  first  five  timesteps.  After  that,  a  rapid  attenuation  of  short  wavelength 
errors  {k^  large)  occurred. 

At  iteration  79,  (Figure  2b),  the  error  is  reset  to  1,  and  the  “fictitious  timestep”  from  the  78th  iteration 
of  previous  problem  is  used  to  restart  the  iterations.  In  this  case,  it  is  the  error  at  long  wavelengths  which 
are  removed  fastest  at  first,  and  some  amount  of  short  wavelength  error  remains  even  after  40  iterations. 
Evidently,  the  “fictitious  timestep”  restarted  at  such  a  value  as  to  remove  the  long  wavelengths  most  quickly. 
The  initial  convergence  is  somewhat  better,  although  the  error  is  not  reduced  as  much  after  40  iterations. 
In  this  example,  changing  the  initial  timestep  improved  the  convergence  in  early  iterations,  but  converged 
more  slowly  at  the  40th  iteration. 
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(a)  40  iterations  on  a  delta  function 


\f(r)\ 


(b)  40  iterations  on  a  delta  function  after  reset 


Iteration 


lE-08 


120  1.0 


Figure  2:  The  effect  of  40  iterations  of  DADI  on  the  error  f{k).  fe'  =  /cAx/tt  is  a  normalized  wavenumber, 
(a)  uses  a  non-optimal  choice  of  initial  “fictitious  timestep”,  which  is  improved  in  (b). 

2.1.2  Langdon-Marder  Correction  on  the  Delta  Function  Source 

Next,  the  Langdon-Marder  correction  is  applied  to  the  delta  function  source.  Longer  wavelengths  are  removed 
less  quickly  by  the  Langdon-Marder  correction  than  the  Boris  correction  (comparing  Figs,  2  and  3),  while 
shorter  wavelengths  are  removed  much  more  efficiently  when  the  numerical  parameter  d  is  chosen  to  be  close 
to  the  maximum  value  possible  for  stability  of  the  algorithm. 

Nevertheless,  even  with  the  choice  of  a  large  d,  the  overall  magnitude  of  the  error  does  not  drop  as 
quickly  (compared  to  the  DADI-based  Boris  correction)  per  iteration.  This  is  balanced  by  the  fact  that  each 
Langdon-Marder  iteration  is  much  less  expensive  than  each  DADI-based  Boris  iteration. 

The  Langdon-Marder  correction  monotonically  reduces  the  error  to  zero  (Langdon  shows  that  the  Langdon- 
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Marder  correction  is  equivalent  to  point  Jacobi  iteration,  which  is  known  to  converge  [3]),  while  the  DADI- 
based  Boris  correction  does  not  converge  monotonically. 

The  choice  of  diffusion  parameter  d  also  seems  to  influence  the  wavelength  at  which  fastest  convergence 
occurs.  However,  very  long  wavelength  errors  are  still  attenuated  slowly:  150  iterations  were  needed  in  some 
cases  to  reduce  the  error  at  the  longest  wavelength  by  a  factor  of  2  with  a  diffusion  parameter  of  d  =  dmax- 


(a)  60  Iterations  of  Marder  (d  =  OJdmax) 
\f(k^)\ 


(b)  60  Iterations  of  Marder  (d  =  0.8  dmax) 


Figure  3:  Attenuation  of  error  as  the  Langdon-Marder  correction  is  iterated,  k'  =  kAx/ir  is  a  normalized 
wave  number,  (a)  was  computed  using  d  =  OAdmaxj  where  dmax  is  defined  by  (7).  The  fiat  part  of  the  graph 
at  10“^  is  an  artifact  of  the  graphing  program;  the  actual  value  of  the  error  there  is  smaller  than  10~®. 
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2.2  Single  Wavelength  Test  Cases 

It  is  also  of  interest  to  examine  the  spectral  diffusivity  properties  of  these  methods.  This  can  be  done  by 
adding  a  divergence  error  at  a  single  wavenumber  and  observing  the  spectral  purity  after  repeated  applications 
of  the  method. 

Accordingly,  error  test  cases  of  the  form 

V-D-p  =  sm(^)sm(^)  (11) 

Lx  Ly 

were  examined,  where  n  is  an  integer  and  Lx,  Ly  are  the  dimensions  of  the  system.  This  type  of  error  shows 
up  as  discrete  spikes  in  Fourier  transforms,  since  only  a  single  wavenumber  is  introduced.  The  values  of  n 
examined  were  chosen  so  that  1,  2,  4,  8,  and  16  waves  could  be  observed  in  the  system. 

For  both  the  Langdon-Marder  and  Boris  corrections,  the  spikes  in  fc-space  remained  discrete  and  did  not 
spread  out.  The  error  remained  in  one  wavelength  until  the  error  was  of  the  order  of  machine  precision. 
Even  errors  with  multiple  wavelengths  introduced  at  the  outset  had  this  property,  although  each  wavelength 
of  error  was  removed  at  a  different  rate.  Since  all  observable  wavelengths  are  present  in  Figs.  2  and  3,  one 
can  observe  the  convergence  rates  by  wavelength. 
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2.3  Relativistic  Beam  Test  Case 

The  delta  function  source  test  case  and  the  single-wavelength  source  test  cases  were  somewhat  contrived, 
in  that  a  divergence  error  was  simply  introduced  into  the  system  arbitrarily.  Those  sources  were  studied  to 
provide  insight  into  the  spectral  behavior  of  the  different  corrections. 

It  is  also  desirable  to  examine  cases  in  which  the  error  arises  from  sources  natural  to  simulations:  the 
actual  errors  introduced  by  moving  simulated  particles.  In  this  case,  a  simulated  beam  with  a  very  high 
velocity  (0.967  c)  was  injected  into  the  system  along  the  x  direction.  The  beam  had  a  low  current  (/  =  lA, 
with  a  crossection  of  1  square  meter)  and  a  high  energy  such  that  the  particles  drift  unperturbed  through 
the  system;  i.e,  gEAt/m  <  y/2£jm  where  E  is  the  initial  beam  energy  and  E  is  the  electric  field.  These 
errors  are  different  in  character  than  the  test-cases  previously  shown,  as  is  evident  from  Figures  4  and  5. 

This  example  is  interesting  because  it  can  generate  large  errors  every  timestep  for  a  current  weighting 
which  is  non-charge  conserving.  The  error  is  shown  in  Figure  4.  Notable  is  the  presence  of  long-wavelength 
components  in  the  error,  as  is  evident  from  Figure  5. 

in  Figure  6,  we  show  the  effect  of  applying  the  corrections.  In  this  figure,  two  iterations  of  the  DADI- 
baised  Boris  correction  enjoy  comparable  success  to  10  iterations  of  the  Langdon-Maxder  correction.  However, 
in  the  top  graph,  the  DADI  method  is  much  more  successful  in  reducing  the  error  in  longer  wavelengths 
but  less  successful  in  shorter  wavelengths  than  the  Langdon-Marder  correction.  In  the  bottom  graph,  both 
methods  enjoy  comparable  success. 

Figure  7  shows  an  overall  figure  of  merit  for  the  removal  of  the  error.  This  figure  was  calculated  as 
follows: 

V  f  /  D  •  dS  -  QoeH))  /  E  (12) 

'  cells 

The  Langdon-Marder  and  DADI-based  Boris  correction  performed  comparably  well.  Two  DADI  iterations 
were  used  for  the  Boris  correction,  and  10  iterations  using  d  =  O.Sdmax  for  the  Langdon-Marder  correction. 
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2.4  Beam  with  a  large  perpendicular  magnetic  field  applied 

This  case  is  similar  to  the  previous  case  except  a  magnetic  field  (63  Gauss)  was  applied  in  the  z  direction  to 
turn  the  beam,  and  a  beam  was  injected  from  the  left  of  the  system  as  well  as  from  the  right.  Both  beams 
were  one  Ampere,  and  travelling  with  a  velocity  of  0.967c.  Figure  8  shows  only  the  rightward-moving  beam. 
The  two  beams  pass  through  each  other  with  little  perturbation,  since  gEAt/m  <  where  E  is  the 

electric  field  and  E  is  the  beam  energy. 

The  error  spectra  observed  in  this  case  (not  shown)  were  very  similar  to  the  spectra  observed  for  the 
previous  relativistic  beam  case,  for  the  uncorrected  simulation  as  well  as  for  the  Langdon-Marder  and  Boris 
corrected  simulations.  The  corrections  were  applied  in  the  same  way  as  the  previous  example:  10  passes  of 
Langdon-Maxder  correction  with  d  =  O.&dmax,  and  2  iterations  of  the  DADI-based  Boris  correction.  One 
notable  difference  is  that  there  is  a  stronger  DC  (fc  =  0)  component  than  in  the  previous  case,  which  was 
shown  in  Fig.  5.  The  stronger  DC  component  probably  explains  why  the  DADI  correction  performed 
a  little  better  than  the  Langdon-Marder  correction  in  Figure  9.  As  was  demonstrated  by  Figure  3,  the 
Langdon-Marder  correction  is  not  as  efficient  at  removing  long-wavelength  errors. 


13 


2.5  Thermal  Plasma 


In  this  case  a  simulation  of  a  very  hot  uniform  plasma  (T  =  4700eV,  rii  =  rie  =  is  examined.  At 

this  temperature  the  electron  thermal  speed  is  3  x  lO^m/s.  The  plasma  is  uniform  with  very  heavy  (M/m 
=  40000)  ions  at  T  =  1.87eV. 

Figure  10  shows  the  time  dependence  of  the  magnitude  of  the  divergence  error  for  the  cases  of  (1) 
uncorrected,  (2)  Langdon-Marder-corrected,  and  (3)  DADI/Boris  corrected.  The  same  number  of  iterations 
and  diffusion  parameter  were  used  as  in  the  previous  cases:  10  Langdon-Marder  iterations  at  d  =  O.Sdmai 
and  2  iterations  of  DADI.  The  DADI-based  Boris  correction  was  again  a  little  more  successful  at  removing 
overall  error  in  the  system.  Again,  the  error  spectra  is  similar  to  that  already  given  in  Fig.  5. 
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3  Conclusions 


Tests  on  single-wavelength  errors  showed  that  both  methods  of  removing  errors  are  non-dispersive,  that  is, 
an  error  at  one  wa-velength  will  not  spread  to  other  wavelengths  as  the  corrections  are  applied. 

Divergence  errors  are  not  necessarily  dominated  by  short  wavelength  errors  for  PIC  simulations.  Instead, 
errors  were  observed  over  a  broad  range  of  wavelengths.  This  broad  spectrum  reduces  the  efficacy  of  methods 
such  as  the  Langdon-Marder  method,  which  preferentially  remove  short  wavelength  errors. 

The  magnitude  of  the  diffusion  parameter  d  used  by  the  Langdon-Marder-type  correction  also  had  some 
effect  on  the  convergence  rate  for  different  wavelengths,  although  this  effect  was  not  very  pronounced. 
Essentially,  the  k  of  fastest  attenuation  of  error  varied  somewhat  depending  on  the  value  of  d:  the  closer 
d  was  chosen  to  dmax,  the  shorter  the  wavelength  of  fastest  attenuation.  However,  even  with  d  near  dmax, 
attenuation  of  long  wavelengths  was  slow. 

The  divergence  error  injected  per  timestep  can  be  large  when  the  particles  are  moving  very  quickly,  and 
thus  error  can  be  introduced  even  more  quickly  than  it  can  be  removed  by  some  corrections.  In  these  cases, 
a  non-charge  conserving  current  algorithm  may  not  be  practical. 

The  divergence  error  was  not  observed  to  increase  without  bound,  but  rather  reached  a  final  value 
asymptotically,  on  the  order  of  1/100  to  1/1000. 

The  Boris  correction  based  on  the  Dynamic  ADI  method  proved  to  be  as  efficient  at  removing  errors  in  the 
PIC  simulations  examined  as  the  Langdon-Marder  correction  used.  Five  iterations  of  the  Langdon-Maxder 
correction  is  about  as  expensive  computationally  as  one  iteration  of  the  DADI-based  Boris  correction  as 
implemented  in  XOOPIC.  Unfortunately,  they  both  produce  almost  the  same  amount  of  correction  overall, 
so  one  cannot  choose  the  better  method  based  on  efficiency  in  error  removal. 

The  Langdon-Marder  correction  is  much  more  straighforwaxd  to  implement  than  a  Poisson  solver,  but 
a  Poisson  solver  may  already  exist  in  the  code  for  other  purposes,  e.g.,  an  electrostatic  field  solver.  The 
Langdon-Marder  correction  is  clearly  easier  to  adapt  to  parallel  processing  architectures  due  to  the  local 
nature  of  the  algorithm.  However,  the  Boris  correction  is  simply  a  matrix  inversion,  amenable  to  library 
routines.  Both  methods  will  scale  as  the  number  of  nodes  in  the  mathematical  mesh. 
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A  How  Errors  in  Gauss’s  Law  Arise  in  PIC  Simulations 


To  fit  the  discrete  mathematical  mesh,  we  recast  V  •  D  =  p  •  dS  =  fy  pdV  around  each  grid  node  as 

shown  in  Figure  11. 

A.l  One-dimensional  non-conservation  of  charge 

Consider  this  one  dimensional  example  of  how  linear  weighting  of  current  and  charge  could  lead  to  a  violation 
of  continuity,  which  in  turn  leads  to  a  violation  of  Gauss’s  law.  In  this  example,  a  particle  of  charge  Q  moves 
from  position  pi  to  p2,  in  a  single  timestep  of  length  At. 

Assuming  this  is  the  only  particle  in  motion,  this  would  lead  to  the  accumulation  of  the  following  charges 
at  mesh  points,  using  a  linear  weighting  of  the  charge  to  adjacent  mesh  points: 
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where  pi,  p2  are  initial  and  final  particle  positions,  Xj  are  grid  positions,  Q  is  the  charge  of  the  particle 
moving,  qj  are  the  charges  weighted  to  the  mesh,  t  is  the  time,  and  At  is  the  time  it  took  for  the  particle  to 
move  from  pi  to  p2- 

Linear  weighting  of  the  current  is  done  by  first  finding  the  midpoint  of  the  displacement  of  the  particle, 
and  then  using  that  midpoint  to  weight  the  current  to  the  nodes  (as  in  [1]  pp.  358-364).  Then,  the  current 
at  the  nodes  is  averaged  to  obtain  the  current  at  the  half-nodes. 
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We  sum  the  currents  into  and  out  of  each  node  to  calculate  the  change  in  charge: 


After  performing  some  algebra,  we  have: 
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Equations  (30)  and  (31)  disagree  with  Eqs.  (17)  and  (18).  Since  the  currents  are  the  source 

currents  used  to  advance  Maxwell’s  curl  equations,  V  x  H  =  J  +  and  V  x  E  =  — 5B/9t,  the  resultant 

electric  fields  refiect  a  different  change  in  charge  at  the  nodes,  i.e.,  a  violation  of  continuity  leading  to  a 
violation  of  Gauss’s  law. 


Similar  difficulties  arise  with  weighting  currents  directly  to  the  half-cells  as  in  Figure  13. 


This  weighting  results  in  the  following  expressions: 
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which  again  differ  from  Equations  (17)  and  (18). 

However,  not  all  current  weightings  will  disagree  with  the  linear  weighting  of  charge  to  the  grid.  If  one 
weights  all  the  current  to  the  nearest  half  grid  point[l,  5],  then  the  current  will  agree  with  Equations  (17) 
and  (18).  The  following  equations  show  this:  (refer  to  Figure  13  for  definitions.) 
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(40) 
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Figure  14  illustrates  the  distribution  of  charge  to  the  grid  for  the  different  current  weightings.  The  figure 
shows  the  integrated  net  charge  at  a  node  due  to  current  weighted  from  a  particle  of  constant  velocity.  At 
time  200  the  particle  is  exactly  on  the  node  point  being  observed.  The  charge-conserving  current  weighting 
produces  charges  which  agree  with  linear  weighting  of  charge.  Linear  weighting  to  half  cells  produces  a  net 
charge  which  agrees  with  quadratic  weighting  of  charge  in  one  dimension,  while  the  linear  weighting  to  nodes 
followed  by  averaging  to  half  cells  does  not  correspond  to  any  straightforward  weighting  of  charge.  Figure  15 
shows  the  current  weighted  to‘  a  half-cell,  which  is  used  in  Maxwell’s  curl  equations  to  advance  the  fields.  At 
time  150  the  particle  is  located  at  the  half  cell  being  observed.  The  area  under  each  curve,  f  Idt,  is  identical 
for  each  weighting  technique.  Figure  15  may  also  be  regarded  as  a  diagram  of  the  shape  of  a  particle  for 
current  weighting  purposes. 

The  disagreement  between  the  charge  given  by  the  linear  and  charge-conserving  weightings  in  Figure  14 
is  the  error,  and  a  source  of  the  violation  of  Gauss’s  law  in  electromagnetic  particle-in-cell  simulations. 

A. 2  Two-dimensional  non-conservation  of  charge 

In  two  dimensions,  the  errors  in  charge  that  occur  with  the  use  of  a  non-charge-conserving  current  weighting 
are  even  more  severe.  Consider  a  bilinear  current  weighting  to  the  half-nodes.  A  particle  traveling  obliquely 
to  the  mesh  can  leave  behind  a  trail  of  permanent  charges  behind  it  due  to  flaws  in  the  weighting.  These 
residual  charges  sum  to  zero,  but  multipoles  are  left  behind. 

As  an  example,  consider  the  charge  on  node  A  in  Figure  16  due  to  a  particle  traveling  uniformly  along 
the  line  abcde.  The  only  weighted  currents  due  to  this  particle  which  affect  the  charge  at  node  A  are  lx  a 
and  lyA‘  The  mesh  is  orthogonal  and  square,  with  Ax  =  Ay.  The  particle’s  velocity  Vx  in  x  and  Vx/2  in  y. 

Only  when  the  particle  is  traveling  from  point  b  to  point  c  will  any  current  weight  to  IxA>  (The  x 
component  of  the  current  due  to  the  particle  is  partially  weighted  to  the  half-node  lx  a-)  The  net  charge 
pushed  into  node. A  due  to  Ixa  is: 


,  qvx  X  /  Ax .  .  1 , ,  1 . 


(42) 


where  qvx/Ax  is  the  current  ,  Ax/2vx  is  the  time  to  traverse  segment  fee,  1/4  is  the  weighting  in  the  x 
direction  to  Ix^,  and  1/8  is  the  weighting  in  the  y  direction. 


Similarly,  weighting  to  lyA  from  segment  ce  yields  a  charge: 


(43) 
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Even  considering  that  =  2vy,  there  is  a  net  charge  left  on  node  A  of  -3q/64. 

Other  nodes  around  the  moving  particle  will  exhibit  similar  errors,  and  another  particle  following  in 
the  same  track  will  double  all  the  errors,  leading  to  charge  accumulation  which  could  severely  affect  the 
simulation. 
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Charge  Density 


Divergence  Error 


Figure  4:  Charge  density  and  divergence  error  (in  C/m^)  in  a  relativistic  beam  simulation,  with  no  corrections 
applied.  Note  that  the  magnitude  of  the  error  is  about  1/4  of  the  magnitude  of  the  charge  density,  i.e., 
Sp/p  ~  1/4.  Both  graphs  axe  noisy  because  of  the  small  number  of  particles  used. 


Uncorrected  error  spectrum 


lError(Kx)!  averaged  over  y 


iError(Ky)l  averaged  over  x 


Figure  5:  Error  spectra  in  a  simulation  with  no  corrections  applied.  Note  that  not  all  the  error  is  in  short 
wavelengths.  To  obtain  better  statistics,  the  error  was  averaged  in  the  direction  opposite  the  FFT. 
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Corrections  applied  to  the  Relativistic  Beam 


IError(Kx)(  averaged  over  y 


Kx 

IEtTor(Ky)l  averaged  over  x 


Ky 


-  Marder 

- Boris-DADl 


-  Marder 

- Boris-DADl 


Figure  6:  Error  spectra  after  corrections  are  applied.  10  Langdon-Marder  iterations  were  used,  with  d 
O-Sdynaxj  while  2  DADI  iterations  were  used.  The  error  was  reduced  by  a  factor  of  more  than  100. 
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Normalized  Divergence  Error 


Simulation  time 


Harder 
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Figure  7:  Normalized  total  error  in  the  simulation. 
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Figure  8:  A  beam  curving  in  a  magnetic  field.  Not  shown  is  the  leftward  moving  beam  entering  from  the 
right  hand  side  and  curving  upward.  Such  a  beam  is  actually  present  in  the  simulation. 
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Figure  13:  Current  weighting  to  half-cells 


Figure  14:  Net  charge  at  a  mesh  point  due  to  weighted  current  from  a  particle  passing  at  a  uniform  velocity. 
The  particle  is  on  the  mesh  point  at  time  200.  The  two  linear  weightings  distribute  charge  differently  than 
the  charge-conserving  weighting,  leading  to  error.  See  the  text  for  details. 
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Figure  15:  Current  at  a  half-cell  due  to  a  particle  passing  at  a  uniform  velocity.  At  time  150  the  particle  is 
located  at  the  half  cell  being  observed. 


Figure  16:  A  particle  traveling  along  abcde  will  leave  a  permanent  charge  at  node  A  if  a  bilinear  current 
weighting  to  half-cells  is  used. 
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2.  Total  funding  of  the  Parent  Agreement  and  the  number  of  full-time  equivalent  graduate  students 
(FTEGS)  supported  by  the  Parent  Agreement  during  the  1 2-month  period  prior  to  the  AASERT 
award  date. 

a.  Funding:  $  125,000 _ 

b.  Number  FTEGS:  One 


3.  Total  funding  of  the  Parent  Agreement  and  the  number  of  FTEGS  supported  by  the  Parent 
Agreement  during  the  current  12-month  reporting  period. 

a.  Funding:  $  505,000 _ 

b.  Number  FTEGS  _ none 


4.  Total  AASERT  funding  and  the  number  of  FTEGS  and  undergraduate  students  (UGS)  supported 
by  AASERT  funds  during  the  current  12-month  reporting  period. 

a.  Funding:  $  137,750 _ 

b.  Number  FTEGS:  _ one 


c.  Number  UGS:  n/a 


VERIFICATION  STATEMENT:  I  hereby  verify  that  all  students  supported  by  the  AASERT  award  are 
U.S.  Citizens. 


Principal  Investigator 

University  of  California,  Berlceley 
F49620-94-1-0387 


7/31/97 

Date 
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